/* 
 * File:   siliconlat.hpp
 * Author: mikdebian
 *
 * Created on July 10, 2013, 3:59 AM
 * 
 * Fill in DiamondCubic lattice
 */

#ifndef DIAMONDCUBICLAT_HPP
#define	DIAMONDCUBICLAT_HPP

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "lattice/atomtypes.hpp"
#include "lattice/lattice.hpp"

#define SILICON_ATOMS_PUC 18

#define SILICON_INTERNAL 1

class DiamondCubic {
public:

    /**
     * Constructor
     * @param cellSize number of unit cells in a cube
     * @param scut cutoff radius
     * @param latconst Lattice constant
     */
    DiamondCubic(int cellSize, double scut, double latconst) {
        mCellSize = cellSize;
        mNUnitCells = 1;
        mA = latconst;
        mSCut = scut;
    }

    void fillLattice(Lattice& lat) {
        int sft = 0;
        lat.mNCell = mNUnitCells * SILICON_ATOMS_PUC;
        lat.mLC = mA;
        double xo[3] = {0, 0, 0};
        fillUnitCell(xo, sft, lat, true);
        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                for (int k = 0; k < 3; k++) {
                    if (!((i == 1) && (j == 1) && (k == 1))) {
                        double x[3], str[3];
                        str[0] = (-1 + i) * mA;
                        str[1] = (-1 + j) * mA;
                        str[2] = (-1 + k) * mA;
                        x[0] = xo[0] + str[0];
                        x[1] = xo[1] + str[1];
                        x[2] = xo[2] + str[2];
                        //printf("<%d, %d, %d> stride = %lf, %lf, %lf\n", i, j, k, str[0], str[1], str[2]);
                        fillUnitCell(x, sft, lat, false);
                    }
                }
            }
            lat.mNA = sft;
        }
    }

private:

    void fillUnitCell(double* xo, int &sft, Lattice &lat, bool prime) {
        double y[3];

        /* Corners */
        y[0] = 0, y[1] = 0, y[2] = 0;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.125;

        y[0] = mA, y[1] = 0, y[2] = 0;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.125;

        y[0] = 0, y[1] = mA, y[2] = 0;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.125;

        y[0] = 0, y[1] = 0, y[2] = mA;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.125;

        y[0] = 0, y[1] = mA, y[2] = mA;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.125;

        y[0] = mA, y[1] = 0, y[2] = mA;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.125;

        y[0] = mA, y[1] = mA, y[2] = 0;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.125;

        y[0] = mA, y[1] = mA, y[2] = mA;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.125;

        /* Facet centers*/
        y[0] = 0, y[1] = 0.5 * mA, y[2] = 0.5 * mA;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.5;

        y[0] = 0.5 * mA, y[1] = 0, y[2] = 0.5 * mA;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.5;

        y[0] = 0.5 * mA, y[1] = 0.5 * mA, y[2] = 0;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.5;

        y[0] = mA, y[1] = 0.5 * mA, y[2] = 0.5 * mA;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.5;

        y[0] = 0.5 * mA, y[1] = mA, y[2] = 0.5 * mA;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.5;

        y[0] = 0.5 * mA, y[1] = 0.5 * mA, y[2] = mA;
        pushAtom(xo, y, sft, lat, prime);
        if (prime)
            lat.mFractions[sft] = 0.5;

        /* Internals */
        y[0] = 0.25 * mA, y[1] = 0.25 * mA, y[2] = 0.25 * mA;
        pushAtom(xo, y, sft, lat, prime, LATTICE_INTERNAL_ATOM);
        if (prime)
            lat.mFractions[sft] = 1;

        y[0] = 0.25 * mA, y[1] = 0.75 * mA, y[2] = 0.75 * mA;
        pushAtom(xo, y, sft, lat, prime, LATTICE_INTERNAL_ATOM);
        if (prime)
            lat.mFractions[sft] = 1;

        y[0] = 0.75 * mA, y[1] = 0.25 * mA, y[2] = 0.75 * mA;
        pushAtom(xo, y, sft, lat, prime, LATTICE_INTERNAL_ATOM);
        if (prime)
            lat.mFractions[sft] = 1;

        y[0] = 0.75 * mA, y[1] = 0.75 * mA, y[2] = 0.25 * mA;
        pushAtom(xo, y, sft, lat, prime, LATTICE_INTERNAL_ATOM);
        if (prime)
            lat.mFractions[sft] = 1;
    }

    void pushAtom(double* xo, double* xn, int &sft, Lattice& lat, bool prime, unsigned int attribute = 0) {
        bool useful = false;
        bool same = false;
        double x[3];
        x[0] = xn[0] + xo[0];
        x[1] = xn[1] + xo[1];
        x[2] = xn[2] + xo[2];
        int limit = (sft < lat.mNCell) ? sft : lat.mNCell;
        for (int i = 0; i < limit; i++) {
            double r = dist(x, lat.mCoor[i]);
            if (r == 0) {
                same = true;
                break;
            }
        }
        if (!same) {
            lat.mCoor[sft][0] = x[0];
            lat.mCoor[sft][1] = x[1];
            lat.mCoor[sft][2] = x[2];
            lat.mTypes[sft] = ATOM_Si;
            lat.mAttributes[sft] = attribute;
            if (prime)
                lat.mNei[sft][0] = 1;
            for (int i = 0; i < limit; i++) {
                double r = dist(x, lat.mCoor[i]);
                //printf("<%d> [%lf, %lf, %lf] - [%lf, %lf, %lf], r = %lf\n",i, x[0], x[1], x[2], lat.mCoor[i][0], lat.mCoor[i][1], lat.mCoor[i][2], r);
                if (r <= mSCut) {
                    useful = true;
                    lat.mNei[i][lat.mNei[i][0]] = sft;
                    lat.mNei[i][0]++;
                    if (prime) {
                        lat.mNei[sft][lat.mNei[sft][0]] = i;
                        lat.mNei[sft][0]++;
                    }
                }
            }
            if (prime || useful) {
                sft++;
            }
        }
    }

    double dist(double* x, double* y) {
        double r = (x[0] - y[0])*(x[0] - y[0]) + (x[1] - y[1])*(x[1] - y[1]) + (x[2] - y[2])*(x[2] - y[2]);
        return sqrt(r);
    }
    int mCellSize;
    int mNUnitCells;
    double mA;
    double mSCut;
};
#endif	/* SILICONLAT_HPP */

